NSF-KITP-03-58 



First-order quasilinear canonical representation of the characteristic formulation of 

the Einstein equations 

Roberto Gomefl 

Pittsburgh Supercomputing Center, 4400 Fifth Ave, Pittsburgh, PA 15213 and 
Department of Physics and Astronomy, University of Pittsburgh, Pittsburgh, PA 15260 

Simonetta Frittellfl 
Department of Physics, Duquesne University, Pittsburgh, PA 15282 and 
Department of Physics and Astronomy, University of Pittsburgh, Pittsburgh, PA 15260 

We prescribe a choice of 18 variables in all that casts the equations of the fully nonlinear charac- 
teristic formulation of general relativity in first-order quasi-linear canonical form. At the analytical 
level, a formulation of this type allows us to make concrete statements about existence of solutions. 
In addition, it offers concrete advantages for numerical applications as it now becomes possible to 
incorporate advanced numerical techniques for first order systems, which had thus far not been ap- 
plicable to the characteristic problem of the Einstein equations, as well as in providing a framework 
for a unified treatment of the vacuum and matter problems. This is of relevance to the accurate 
simulation of gravitational waves emitted in astrophysical scenarios such as stellar core collapse. 
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I. INTRODUCTION 

The characteristic formulation of general relativity due 
to Bondi and Sachs 0, Q, in particular the approach 
based on a null slicing of spacetime with a transverse 
timelike data surface as proposed by Tamburino and 
Winicour has been used successfully for many ap- 
plications in numerical relativity. It provides a natural 
framework for the computation of gravitational radia- 
tion signals from isolated astrophysical sources by purely 
characteristic evolution 4\ , and as a means of extracting 
gravitational radiation information from 3 + 1 simula- 
tions 

It has been used to achieve long-term stable numerical 
evolutions of generic single black hole spacetimes 0, Q , 
it has also been used to compute the behavior of matter 
fields around black holes and is ideally suited to 

the study of gravitational radiation of black hole-neutron 
star mergers, the prime candidates for detection by ad- 
vanced gravitational wave interferometers |ll| . A modi- 
fication proposed recently incorporates into the for- 
mulation a partial treatment of matter fields with the 
specific goal of modeling the capture of a massive ob- 
ject such as a neutron star by a galactic-size black hole, 
up to the point where tidal disruptions become impor- 
tant. Events of this type are the predominant sources of 
gravitational radiation which are expected to fall in the 
frequency band of LISA. 

Additionally, the formulation provides a unique ap- 
proach to the post-merger regime of binary black hole 
coalescence, starting from the gravitational radiation 
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emitted during a white hole fission as illustrated 

in [T3. Hfij in the close-limit of the binary black hole 
merger. It has also been used to study the nonlin- 
ear generation of waveforms in single black-hole space- 
times 0, 0| , where it has shown the potential to gen- 
erate a catalog of waveforms, also of interest for data 
analysis of space-based gravitational interferometers. 

One would expect the stability properties exhibited 
by numerical representations of the characteristic prob- 
lem 0, 0, IE S t° reflect underlying stability prop- 
erties at the analytical level, perhaps along the lines of 
[T^.l20| . In this regard, in [2lJ the linearized equations are 
cast into a canonical first-order form that is suitable for 
the use of Duff's theorem of existence of solutions |2^ . 
However, the choice of variables of [2l| does not allow 
for an extension of the result to the non-linear case in 
any obvious manner. It is desirable to have a quasilinear 
formulation of the characteristic equations 0, 181 in first- 
order form that is to the characteristic problem what a 
first-order formulation is to a Cauchy problem. Such 
a formulation could in principle be used as the starting 
point to approach relevant issues of stability by means of 
energy estimates. 

Of particular importance to the numerical implemen- 
tation of the characteristic approach is that by writing 
the system of equations in first-order quasilinear form, 
conservative schemes and Godunov-type shock capturing 
methods |23l l24j commonly employed to solve non-linear 
hyperbolic systems of conservation laws (e.g. the Euler 
equations) can be brought to bear on the gravitational 
field equations [25j. High resolution shock capturing 
methods have been successfully used in the characteristic 
formulation for the evolution of matter fields 0, , but 
they have not been applied to the equations for the grav- 
itational fields themselves, neither in the 3-dimcnsional 
case, with the equations in the form presented in Ref. 
or in Ref. Q, nor, to the best of our knowledge, to the 
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equations for the axisymmetric case |26(. The approach 
developed by Pons et. al [53] extends the use of special 
relativistic Riemann solvers, developed for special rela- 
tivistic flows |2^, through local coordinate transforma- 
tions, thereby making them applicable to general rela- 
tivistic problems as well. An obvious impediment to a 
unified treatment, as proposed in Pons et. al. ^2%, has 
been the lack of a first order representation of the char- 
acteristic formulation. 

This leads to a subtle point deserving of clarification: 
true shocks do not form in the gravitational field, thus 
strictly speaking, shock capturing methods are necessary 
for the matter evolution equations only. In the presence 
of matter however, given that the matter fields act as 
sources for the metric fields, shocks and steep gradients 
in the matter sources lead invariably to the appearance 
of short-scale spatial features in the gravitational fields. 
The centered finite difference schemes used in the char- 
acteristic codes 0, 0, 0, |2(| are not capable of resolving 
these short scales, and invariably generate high frequency 
noise. In a perverse twist, the more successful a shock 
capturing method is at resolving short-scale matter fea- 
tures, the more it compounds the problem of integrating 
the metric fields across these features by means of stan- 
dard centered difference schemes. The effect is a classical 
illustration of the error one finds in advecting a pulse us- 
ing a non-shock capturing centered difference method, 
e.g. Lax-Wendroff, where spurious oscillations appear 
on the trailing edge of the numerical solution. With a 
convergent scheme, these unwanted oscillations decrease 
with increasing resolution and vanish in the continuum 
limit, but are present in all simulations with practical 
grid sizes. Attempts to filter out the noise lead to a 
spreading of the pulse. In the problem at hand, it can 
lead to an unacceptable trade-off between the accuracy 
of the matter fields evolution vs. that of the metric 
fields. This effect is mentioned in recent work by Siebel 
et.al. '291, who perform characteristic evolution in ax- 
isymmetric stellar core collapse. In their work, the mat- 
ter fields are solved via shock capturing methods while 
the metric evolution is treated by centered, second order 
accurate finite difference schemes |26| , with modifications 
along the lines of . This is a difficult numerical task: 
in the core collapse scenario, a strong shock wave forms 
after bounce, where all matter fields are discontinuous. 
The dynamics of the collapse are correctly solved and 
these discontinuities are accurately tracked [29| . How- 
ever, the discontinuities in the matter field lead to un- 
avoidable discontinuities in the first derivatives of the 
metric fields. In the traditional scheme, second deriva- 
tives of the metric are need to compute the gravitational 
radiation, and this is where numerical noise is generated, 
the main difficulty pointed out in |29j. A form of the 
equations without second order derivatives would pro- 
vide a solution to this problem, and thus be particularly 
relevant to astrophysical problems with shocks. In the 
absence of such a formulation, filtering of the unwanted 
noise on a fixed grid has been shown not to provide suf- 



ficient accuracy, and the conclusion drawn in 29] is that 
unless the short scale gravitational field features are ad- 
equately resolved and the resulting spurious oscillations 
eliminated, it is difficult to extract more accurate gravi- 
tational signals, only the main features of the signal being 
obtained to good accuracy. 

It is generally accepted that it will be necessary to 
incorporate adaptive mesh refinement techniques in or- 
der to achieve well-resolved simulations with matter 
sources 113 • Preliminary steps in these direction 
have been taken, with the implementation of an adap- 
tively refined code for the model problem of Einstein- 
Klein-Gordon fields in spherical symmetry 30] . This task 
can also be more easily addressed when both the gravita- 
tional field equations and the matter equations are writ- 
ten in first order differential form. We should also point 
out another difference between the matter and gravita- 
tional field equations. Since the matter evolution equa- 
tions express the conservation of the stress-energy ten- 
sor, they can be easily expressed in conservative form 
The gravitational field evolution equations, on the other 
hand, need not be put in conservative form, as they can 
be treated with volume preserving methods which ap- 
ply equally to first order quasilinear systems and are also 
available |3I| . 

Here we introduce a set of variables and auxiliary equa- 
tions that casts the nonlinear gravitational field equa- 
tions into a quasilinear, first-order representation of the 
Einstein equations in the null cone formalism that takes 
Duff's canonical form and provides a bridge to poten- 
tial adaptations of Cauchy methods to the characteris- 
tic problem for the Einstein equations. The Tamburino- 
Winicour version of the Bondi-Sachs characteristic prob- 
lem reduced to first-order in angular derivatives is re- 
viewed in Section HTIand is taken as the starting point for 
the remainder of the work. In Section ITTT1 the first-order 
quasilinear canonical characteristic form of the equations 
is derived. We offer concluding remarks in Section llVI 



II. THE CHARACTERISTIC PROBLEM OF 
THE EINSTEIN EQUATIONS 

As in Refs. HHEl, we use coordinates based upon 
a family of outgoing null hypersurfaces. We let u label 
these hypersurfaces, x A (A — 2, 3) label the null rays 
and r be a surface area coordinate, such that in the x a 
(it, r, x A ) coordinates the metric takes the Bondi-Sachs 
form 0,0 



ds 2 



r 2 h AB U A U B ) du 2 - 2e 2l3 dudr 

r 



2r 2 h AB U B dudx A + r 2 h AB dx A dx B , (I) 



where h A B is conformal to the metric of the sections 
of fixed value of r on the null slice, and det(h A B) = 
det(q A B), with q A s a unit sphere metric. We define 
the inverse by h AB hsc — $&• By representing ten- 
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sors in terms of spin- weighted variables the confor- 
mal metric Hab is completely encoded in the complex 
function J = \\iABq A q B , where q A is a complex dyad 
such that q A q~A = —2 and q A qA = (an overbar de- 
notes complex conjugation). The remaining dyad com- 
ponent of the conformal metric, given by the real func- 
tion K — \h,ABq A q B is determined by K 2 = 1 + J J as 
a consequence of the determinant condition. Addition- 
ally, we define U = V \q A ' ■ Angular derivatives of tensor 
components are in turn expressed in terms of S and <3 
operators 32]. 

The equations for the characteristic (or null cone) for- 
mulation follow from projections of the Ricci tensor nor- 
mal and tangent to the null slices 0] • The resulting main 
equations arrange into a hierarchy, splitting into a set of 
hypersurface equations, which involve only derivatives on 
the null cone, and evolution equations involving deriva- 
tives with respect to the retarded time u. In particu- 
lar, R rr — provides an equation for /3 r in terms of J, 
while R r Aq A = gives U, rr in terms of J and (3, and 
the trace Rab^ ab = yields V r in terms of J, (3 and 
U . Finally, RABq A q B supplies the evolution equation 
for J. The remaining four components of the Ricci ten- 
sor vanish as a consequence of these six in the following 
sense. By virtue of the Bianchi identities, the compo- 
nent R ur = is trivially satisfied wherever the main 
equations are satisfied, whereas the remaining compo- 
nents R uu = and R u Aq A — are propagated radially 
on the null slices if they hold on a surface r = r$. Thus 
Ruu = and R u Aq A = (the supplementary conditions) 
can be viewed as constraints on the data at r = r$. In 
the following, we ignore these constraints. We will also 
ignore matter source terms, although including them is 
straightforward (see Q), in the interest of keeping the 
presentation concise. 

For the present derivation, we find it convenient to 
start from a relatively recent representation of the char- 
acteristic formulation 0] which casts the system into 
first-order form in the angular derivatives, and mixed 
first-second-order form in the radial derivatives -as op- 
posed to the standard, mixed-order form 0, . We de- 
part however from in our choice of fundamental vari- 
ables in order to facilitate the presentation in Sec. IIIII In 
this partially reduced form, the complete system of main 
equations of the characteristic formulation consists of a 
complex evolution equation for the conformal metric 

2 (rJ), ur - ((! + r W) (rJ) >r ) r = V + J H + JP U , (2) 

namely Eq. (16) of Ref. Jig, and the hierarchy of hyper- 
surface equations and auxiliary definitions as follows: 



{r 2 Q),r 



2r 2 B r - ArB + r 2 



K(k jr + v /r ) 



+ vj^ r + Jfl. r 



uK, 



r 2 U, 



+ — r [P(j r - J 2 J r ) + fl(J >r - J 2 J, r )] (7) 

(8) 



(r 2 W), r 



2K 2 

= e 2fi (KQ-JQ) 
1 



2 e 2f3 TZ - 1 + rW + rW + — (SC/ r 



+ e 
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(— K (pB + BB) 



-[J (dB + B 2 ) + J(BB + B 2 )} 

+~[(v-k)B + (p-k)B] 
r 4 r- 

e' 2P — U, r (KU.r + JU, r 

8 L 

+ U r (KU r + JU r 



(9) 



These are Eqs.(21)-(26) of [l§| with the identification 
\i = 8 J. In Eq. the left-hand side is a characteristic 
representation of a wave operator of second differential 
order in r and u, and we have split the right hand-side 
into three terms, according to usage. The symbol P u sin- 
gles out the only term where a retarded time derivative 
appears in any right-hand side, namely J, u : 

Pu = ^ [J,u (JrK - JK, r ) + lu (JrK - JK, r )] (10) 

The occurrence of this retarded time derivative in first 
order will be found to be critical to the developments 
that follow in Section ITTT1 The symbol Jh stands for the 
following collection of terms: 

e 213 ( 

,1 H = — ( - K{fiB + 2kB - vB) + B{Jfi + Jv) 
+J(Bk - Bk) +J[- 2K(dB + BB) 
+J(BB + B 2 ) + J(m + B 2 )] 



-20 



{KU, r + JU. r ) 

[(KU : r + JU, r )U r + {KU.r + JU. r )U., 



2 
J 

~2 

~(rU r + 2U)-^(rU, r + 2U) 

+~(r§l7, r + 2W) - ~(rS[/ r + 2W) 
+(1 - K)(rdU r + 2W) - ^J, r (SU + BU) 





= 9J r 


(3) 




V.r 




M,r 


= 9J, r 


(4) 


+ r(J r K- 


Ar 


= 7 -(j r J r -K 2 r ) 


(5) 




B r 


= S#r 


(6) 


+JBU- 



(11) 
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Any remaining terms are collected into the symbol T> for 
notational convenience: 

2<= 2 < s 

V = -K(BrU r +2QU) + (dB + B 2 )-(rW) r J. (12) 

r 

We also have 

n = 2K+ 1 -[3{ v -k) + 3{D-k)} +i l(| M | 2 -H 2 ), 

(13) 

and use the symbol k throughout for 

(») 



The symbol is simply a renaming of the original met- 
ric function V that is regular at r = oo, being defined 



by V 



~W. The symbol Q is a first-order vari- 



able encoding the radial derivative of U, and is defined 
by Eq. (JSJ, which acts as a hypersurface equation for U. 
In addition, z/, /i and B are first-order variables of spin 
weights 1 , 3 and 1 , respectively, used to reduce the differ- 
ential order of the angular derivatives appearing in the 
original characteristic equations, and are defined by 



v = 8 J, fi = 9 J, B = 



(15) 



With the definitions {T1J-{T5j), Eq. {TTJ follows from 
Eq. (25) of 0. Equations ©-© as they stand con- 
tain no second-order derivatives in the retarded time or 
the angular coordinates. Additionally, all appearances of 
U, r and W, r in the right-hand sides represent angular 
derivatives and undifferentiated terms by virtue of Q 
and @. The system is still of second differential order 
overall, because of the presence of second-order deriva- 
tives of J, but its value resides in the fact that it exhibits 
remarkable numerical stability properties , raising the 
question of whether a full reduction to proper first order 
may further enhance the numerical stability. 



III. REDUCTION OF THE 
TAMBURINO-WINICOUR SYSTEM TO 
FIRST-ORDER QUASILINEAR FORM 

Our immediate goal is to write the full nonlinear equa- 
tions in a quasi-linear first-order form in the strict sense, 
that is: 



A a (u, r, x A , v)v >a + s(u, r, x A ,v) = 0, 



(16) 



where v represents the set of all dependent variables, the 
index a runs over all spacetime coordinates and where 
the matrices A a and the vector of source terms s depend 
on the coordinates and the undifferentiated variables v. 
Since all remaining second-order terms contain J, r , and 
because the only nonlinear terms in first- derivatives are 
quadratic and contain J, r as a factor - see Eq. (|10|) — . 
this can be accomplished if an appropriate r— derivative 
of the fundamental field J is re-defined as an additional 



fundamental variable, and there is any number of dif- 
ferent acceptable re-definitions, one of which was used 
in [2l| ■ Proceeding along the lines of Ref. 0] , we define 



H = (rJ), r , 



(17) 



which has spin weight 2. No other radial derivatives 
are necessary to convert Eqs. lHJ-lJlJl down to first-order 
quasilinear form. With this definition, the left-hand side 

of Eq. © becomes 2H, U - ((1 + rW)H^ , r . In the pro- 
cess, however, the term J lU in the right-hand side of 
Eq. (0 is promoted to the principal symbol of the sys- 
tem, with the consequence that the evolution equation in- 
volves the retarded time derivatives of two complex vari- 
ables (H and J), instead of just one. It is unclear at this 
point whether Eq. would determine the evolution of 
H or of J. (The difficulty does not arise if one linearizes 
the equations before reducing to first-order form, as was 
done in |21|.) In fact, we know of no solution-generating 
process for the system at this point. So we proceed as 
follows. 

In order to avoid the occurrence of the retarded-time 
derivative of J in the right-hand side of Eq. (J2J we define 
it as an additional fundamental variable: 



F = Ju, 



(18) 



which has spin weight 2. This is at first counter-intuitive: 
the equations are already first order in d u , so defining the 
u — derivative as a new variable might not appear neces- 
sary, or even consistent. However, in the following we 
show that by defining this additional variable, the char- 
acteristic equations take the canonical hierarchical form 
needed for the existence of a solution from characteris- 
tic data [l^. With the definitions (|T7J> - ljT%|> . the original 
evolution equation, Eq. @ , is interpreted as a wave equa- 
tion for H, shown below as Eq. | |19j| ). From (|18fl we have 
(rF)^ r = (rJ), ur , which yields a hypersurface equation 
for F, namely Eq. I|19if) below. With this we can finally 
write the system in the form 





r 


(19a) 




= -(3ff-M), 
r 


(19b) 




r 


(19c) 


P,r 




(19d) 


8rB r 


= r\i, r {H - J) + rD :T (H - J) 





1 

K 



J(H - J) + J(H - J)] rfc r (19e) 
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(r 2 Q), r = 2r 2 B, 



4rB 



— K{k. r + I/ r ) 



+vj r + J jl t r + fiT r + Jk r — J r k 



r 2 [/, 



(r 2 V^), r 



2K 2 
e 2/3 (KQ ~ JQ) 
1 



1 



rW + r5£7 
i(S [e 2p (KQ - JQ)] +8 
AT (8B + -B-B) 

J(g J B + J B 2 ) + j (8s -i 



2,1 



2(rF),r 



+-[(1/ + 
e 2/3 

— [Q(#Q - JQ) + Q(J4TQ - JQ)] 

o 

(l + rWOfr) +V + J H + JP U , 



J 2 J,r)] (19f) 

(19g) 

[e 2/S (ifQ- JQ)]) 
S 2 )] 



(19h) 
(19i) 



for the hypersurface equations and 

2H U - [(1 + rW)H] >r =V + J H + JP U , (19j) 
for the evolution equation, where 
P u = F(H-J)+F(H-J) 

1 1 + 1 "' ^ [{H - J) J + J(H - J)] (20) 



2K 2 



Eqs. H19b|) . i|19c|) . (|19e|) and (|19i|) arise from taking an 
r— derivative of l|15f) and commuting the derivatives in 
the right-hand side, as usual. All the radial derivatives 
indicated in the right-hand side of Eqs. (|19el) -( [T^j| can 

be substituted, in turn, by quantities computed previ- 
ously in the hierarchy, i.e. the right-hand sides of the 
equations can be expressed purely in terms of the un- 
differentiated fundamental variables and their angular 
derivatives. The substitutions have been left indicated 
for the sake of brevity, noting only that the derivatives 
of k (which is not part of the hierarchy) are given by 



2K 



kK., 



' jk = 2K + JdD + ^ 



K 



Equations i|19a[l - i|19i|) can be viewed as propagation equa- 
tions along the radially outgoing null geodesies, with 
Eq. d 1 9 j D advancing the radial derivative of the spheri- 
cal metric function J in time. 

In this first-order formulation of the null cone ap- 
proach, Eqs. H19af> - |T§j| ), the boundary data at r = r 
consists of the values of J, (3, Q, U and W, with the val- 
ues of /i, is, B and F following from the boundary value 



of J as per Eqs. I|15(l and (|17|l . The consistency condi- 
tions, imposed at r = r , are propagated to the interior 
by Eqs. (fl9bjl . (jT3c|) . l(T^e|l and $T§$. The initial data 
for the system (|19aH -( fEJj| at u = uo are the values of 
H{r,x A ), representing the shear of the conformal metric 
of the spheres, given on the entire initial hypersurface. 
The conformal metric function J on the initial hypersur- 
face follows by integration of H as per Eq. I|19a|) . with 
the integration constant provided by the value of J at 
the boundary. Eqs. (|19b|) through (|19i|l in turn provide 
initial values for /i, v, k, f3, B, Q, U , W and F, while 
Eq. | |19j| l propagates H forward in retarded time. At 
this point, the process can be repeated, and the entire 
space-time exterior to the time-like data surface can be 
computed. 

This solution-generating process lies at the basis of 
Duff's theorem of existence and uniqueness of solutions 
to characteristic problems -that is: problems for hyper- 
bolic systems of equations where data is prescribed on a 
characteristic surface. From the analytical point of view, 
as it stands, the system of Eqs. i)19f) has the form 

d u q + Nd r q = Jj 1 (9g, Bq, (5w, Bw, q, w) (22) 
d r w + Md r q = L 2 (Bq,Bq,(5w,Bw,q,w) (23) 

where q = H, w = (J, /i, v, (3, B, U, Q, W, F), and N and 
M are certain matrices of dimension 2x2 and 14 x 2 re- 
spectively, depending on the undifferentiated variables. 
A trivial change of variable F — » F — (1 + rW)H/(2r) 
puts the system of equations l|19a|) -( [T§jl ) into a 18- 
dimensional first-order canonical quasi-linear form as de- 
fined by Duff [23, for 18 variables of which two (H) are 
normal and the remaining 16 (J, fi,v, /?, B,U,Q,W , F) 
are null, and with four complex constraints Ci on the sur- 
face r = tq which are trivially preserved by the solution- 
generating process in the form d r Ci = 0. This means that 
the system (|19|) satisfies the conditions for Duff's theo- 
rem, and therefore existence and uniqueness follow from 
null and normal data in a manner analogous to Cauchy 
problems. This is not a trivial result, as readers should 
note that the same cannot be said if J, u is not defined as 
a fundamental variable. 



IV. CONCLUSION 

The system of equations (!19all -( [T^j| ) casts the 
Tamburino-Winicour version of the Bondi-Sachs charac- 
teristic initial value problem into a standard first-order 
quasilinear form. 

A novel feature of this formulation is the introduction 
of a u— derivative of the 2-sphere metric, 1t,ab,u, as a fun- 
damental variable. The necessity of this step arises only 
in the full non-linear characteristic problem, signaling the 
fact that the linearization and "canonization" operations 
(i.e. the reduction to a hierarchy as per the Bondi Sachs 
construction) do not commute in the case of the charac- 
teristic problem of the Einstein equations. This form of 
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the equations opens the possibility of further studies at 
the analytic level, specifically to look for the existence 
of estimates of the solution in terms of the data on the 
initial characteristic surface and the data on the surface 
of fixed radius r^. 

In a much broader context, we have constructively 
shown here that only 18 variables are needed in or- 
der to achieve a first-order formulation of the Einstein 
equations suitable for numerical generation of solutions. 
And, most remarkably, the resulting nonlinearities are of 
the quasilinear type. By contrast, with regards to the 
Cauchy problem of the Einstein equations, Alekseenko 
and Arnold [3j| show that as many as eight variables are 
actually needed in addition to the six three-metric com- 
ponents and the six extrinsic curvature components in 
order to obtain a full first-order reduction of the ADM 
equations £ 34]. This yields a total of at least 20 vari- 
ables for the first-order Cauchy problem, with the draw- 
back that the resulting non-linearities are genuine (not 
of quasilinear type). In order to remove the genuine non- 
linearities, all 18 space-derivatives of the three-metric 
must be added as fundamental variables, with the re- 
sult that the generic quasilinear first-order reduction of 
the 3+1 Einstein equations requires a number of 30 vari- 
ables. From this perspective, the fact that only 18 vari- 
ables in all are actually sufficient for a first-order quasi- 
linear solution-generating process for the Einstein equa- 
tions from given data is both surprising and intriguing. 



But perhaps the most important feature of the first- 
order formulation Eqs. (|19a[l - |lTjj| ) is its potential for ac- 
curately handling discontinuities in the first derivatives of 
the metric. This is very relevant to simulations of systems 
of astrophysical interest involving shock waves, such as 
stellar core collapse, where discontinuities in the matter 
fields arise and are transmitted to the first derivatives 
of the metric. In such systems, an accurate treatment 
of those discontinuities is essential to ensure the qual- 
ity of the waveforms obtained. The quasilinear form of 
the equations is thus ideal for a unified treatment of the 
gravitational and matter evolution equations, with the 
introduction of more advanced numerical algorithms, in 
particular along the lines of 0, El E3i EH • Work in these 
directions is currently in progress. Results of the appli- 
cation of the system of equations introduced here to the 
numerical characteristic effort will be reported elsewhere. 
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